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^T) I Abstract 

In this paper, we consider numerical homogenization of acoustic wave 
"^ I equations with heterogeneous coefficients, namely, when the bulk modulus 

^-i ' and the density of the medium are only bounded. We show that under 

r~| ■ a Cordes type condition the second order derivatives of the solution with 

j^ i respect to harmonic coordinates are L^ (instead H~^ with respect to Eu- 

clidean coordinates) and the solution itself is in L°°(0,T, iJ^(Q)) (instead 
of L^{0,T,H^{Q)) with respect to Euclidean coordinates). Then, we pro- 
ofs . pose an implicit time stepping method to solve the resulted linear system on 
^ I coarse spatial scales, and present error estimates of the method. It follows 
^^ ■ that by pre-computing the associated harmonic coordinates, it is possible 
^f^ . to numerically homogenize the wave equation without assumptions of scale 
'nJ" I separation or ergodicity. 
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1 Introduction and main results 

Let r2 be a bounded and convex domain of class C^ of M^. Let T > 0. Consider 
the following acoustic wave equation 
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' K-\x)d'fu = dw{p-\x)Vu{x,t)) + g in fix(0,r). 

u{x,t) = for {x,t) edQx {0,T). 

u{x,t) = u{x,0) for {x,t) eQx {t = 0}. 
^dtu{x,t) = ut{x,0) for {x,t) eQ X {t = 0}. 

Write Qt := fi x (0, T) and a := p~^. We assume a is a uniformly elliptic 2x2 
symmetric matrix on Q whose entries are bounded and measurable. There exists 

< ajnin < Omax, SUch that V^ G M^, 1^1 = 1, ttmin < ^^a{x)^ < ttmax, Vx G f2. K is 

a scalar such that -K'min ^ K < K'max- 5' ^ -^^(^2^). 

Equation fll.ip can be used to model wave propagation in heterogeneous media. 
It is important in many applications such as geophysics, seismology, and electro- 
magnetics P, ini EOl [32j . In geophysical and seismic prospecting, K stands for the 
bulk modulus, p the density and u the unknown pressure. The velocity c and 
acoustic impedance a are given by c = \fKjp and o = \JKp. 

Wave propagation in heterogeneous media involves many different spatial 
scales. Even with modern state-of-the-art supercomputers, a direct simulation 
of the highly heterogeneous media is often difficult if not impossible. That is why 
we want to use multiscale methods to solve (11.11) on the coarse spatial scales. More 
precisely, we want to know how to transfer information from fine scales to coarse 
scales, and how to use the information obtained to solve the coarse scale problem 
with much fewer degrees of freedom. We often refer this procedure as numerical 
homogenization or numerical upscaling. 

The idea of using oscillating tests functions in relation to homogenization 
can be backtracked to the work of Murat and Tartar on homogenization and H- 
convergence, in particular we refer to [2S] and [31] (recall also that the framework 
of H-convergence is independent from ergodicity or scale separation assumptions). 
The implementation and practical application of oscillating test functions in fi- 
nite element based numerical homogenization have been called multi-scale finite 
element methods (MsFEM) and have been studied by numerous authors [3], [7], 
|12] . [T3], [20]. On the other hand, numerical schemes have been developed to 
solve the acoustic wave equation with discontinuous coefficients, for example in 
[9] by nonconforming finite element method and in [6] by domain decomposition. 
Recently, numerical homogenization or numerical upscaling methods such as [32] 
are proposed for wave equation with heterogeneous coefficients. 

The finite element method in this paper is closer in spirit to the work of Hou 
and Wu [2D] and Allaire and Brizzi [3] . It is based on a technique first introduced in 
[28] for elliptic equations and extended in [27] to parabolic equations characterised 
by a continuum of scales in space and time. The main difference lies in the fact 
that, instead of solving a local cell problem to get a basis function as in MsFEM 
(Multiscale finite element method) or to calculate effective media property as in 



upscaling method [H] , we use a global change of coordinates. The global change 
of coordinates allows to avoid the so called cell resonance problem and obtain 
a scheme converging in situations where the medium has no separation between 
scales. This makes our method amenable to problems with strongly non-local 
medias, such as high conductivity channels. 

We use a composition rule to construct the finite element space. I. Babuska et 
al. introduced the so called "change of variable" technique [7] in the general setting 
of partition of unity method (PUM) with p- version of finite elements. Through this 
change of variable, the original problem is mapped into a new one which can be 
better approximated. Allaire and Brizzi [3] introduced the composition rule in the 
multi-scale finite element formulation, and have observed that a multi-scale finite 
element method with higher order Lagrange polynomials has a higher accuracy. 

The main difference with parabolic equations [27j lies in the fact that with 
hyperbolic equations, energy is conserved and after homogenization there is no 
hope of recovering the energy (or information) lying in the highest frequencies. 
However when the medium is highly heterogeneous the eigenf unctions associated 
to the highest frequencies are localized, thus energy is mainly transported by the 
lowest frequencies. That is why, when one is only interested in the large scale 
transport of energy it is natural to approximate the solutions of (11. ip by the 
solutions of an homogenized operator. For localization of waves in heterogeneous 
media, we refer to [S [231 [22] • 

This paper is organized as follows. In the next section, we present the formu- 
lation of the mathematical problem and numerical methods, and also show main 
results. In Section 3, we will give the detailed proof and explanation of the results 
in Section 2. In Section 4, we present several numerical examples and conclusions. 

2 Main Results 

In general, the approximation power of finite element method is subject to the best 
approximation for an exact solution with respect to the finite element space. There- 
fore, we require smoothness of the solution to prove convergence theorems. That 
is one of the reasons why standard methods are not applicable for problems with 
heterogeneous media. For example, in (11.11) . we only have u G L°°{0,T,H^{Q)), 
and we can not gain anything if we approximate the solution with usual C° or C^ 
finite element basis. However, as in |23], we can find harmonic coordinates which 
the solution of the wave equation is smoothly dependent on, which is the so-called 
compensation phenomena. 



2.1 Compensation Phenomena 

We will focus on space dimension n = 2. The extension to higher dimension 
is straightforward conditioned on the stability of a. Let F := (^1,^2) be the 
harmonic coordinates satisfying 



divaVF = in n, 
F[x) = X on dil. 



(2.1) 



Let a := ^VFaVF and 



/x^ := esssup^gf^ . . j . (2.2) 

Condition 2.1. a satisfies Cordes type condition if: yUg- < 00 and (Trace[o"])^^ G 
L°°(fi). 

Remark 2.1. If F is a quasiregular mapping, i.e., the dilation quotient (the ratio 
of maximal to minimal singular values of the Jacobi matrix) is bounded, then 
Cordes type condition 12.11 is satisfied [2] . A invertible quasiregular mapping is 
called quasiconformal. In [2] and references therein, invertibility of F is proved 
for a G L°°(f2). Some sufficient conditions for F being quasiconformal were also 
given, for example, det(a) is locally Holder continuous. Unfortunately, a coun- 
terexample with checkerboard structure was proposed, and it can be shown that 
/io- is unbounded at the intersecting point, which is known in mechanics as stress 
concentration. However, we will show that as a solution technique, the numerical 
methods proposed in this paper also works for the cases with stress concentration. 

Let L^(0,T; Hq^H)) be the Sobolev space associated to the norm 

rT 
\\'^\\L^(O.TMi(n)) -^ / r(-'^) mro^"^ (2-3) 



t/ u 

Also, we define the norm of the space L°°{0, T, H'^{Q)) by 

\\v\\L^{o,T,H^n)) = esssupo<t<^( / V" {didjv{x,t)Y dx) \ (2.4) 

We require the right hand side g, initial value u{x, 0) and Ut{x, 0) to be smooth 
enough, which is a reasonable assumption in many applications. For example, we 
can made the following assumptions. 

Assumption 2.1. Assume that the g satisfies dtg G L^(r2r), g G L°°(0,T, L^(i7)), 
initial data m(x, 0) and ^(^(a;, 0) satisfy dtu{x,0) G H^{Q) and Va{x)'Vu{x,0) G 
L^(fi) or equivalently d^u{x,0) G L'^{Q). 



We have the following compensation theorem, 

Theorem 2.1. Suppose that Cordes condition \2. 1\ and Assumption \2. 1\ hold, then 
uoF-^ eL^iO,T,H^{n)) and 

\\uoF~^\\Lo.(o^T,min)) <C(||fi'||L°°(o,T,L2(f7)) + WdtgHmnr) + \\dtu{x,0)\\min) 

+ \\dfu{x,0)\\L2(n))- 

The constant C can he written as 

C = C(n,ll,irmin,^max,amm,amax)/Ua||(Trace[cr])"^||^^^^^. (2.6) 

Remark 2.2. We have gained one more order of integrability in the harmonic co- 
ordinates since in general u G L'^{0,T,H^{Q)). The condition g G L'^{Qt) is 
sufficient to obtain Theorem 12.11 and the following theorems. For the sake of clar- 
ity we have preferred to restrict ourselves to (7 G L°°(0,T, L^(0)). 

2.2 Numerical Homogenization in Space 

Suppose we have a quasi-uniform mesh. Let X^ be a finite dimensional subspace 
of Hq{Q) n W^'°°{Q) with the following approximation properties: There exists a 
constant Cx such that 

• Interpolation property, i.e., for all / G H^{Q) fl Hq{Q) 

mi \\f-v\\Hi^^^<Cxh\\f\\HHn). (2.7) 

• Inverse Sobolev inequality, i.e., for all v G X^, 

\\v\\m{n) < Cx/i^l^^||/fi(f7), (2.8) 

and 

||^^||//i(n) < Cxh'^\\v\\L2(n)- (2.9) 

These properties are known to be satisfied when X^ is a C^ finite element space. 
One possibility is to use weighted extended B-splines (WEB) method developed 
by K. Hollig in [181 [IH], these elements are in general C^-continuous. They are 
obtained from tensor products of one dimensional B-spline elements. The homoge- 
neous Dirichlet boundary condition is satisfied by multiplying the basis functions 
with a smooth weight function uj which satisfies cj = at the boundary . 

Write the solution space V^ as 

V'':= {ifoF{x) : (/p G X'^}. (2.10) 



Remark 2.3. We prove all the following theoretical results by using exact F, how- 
ever, in the numerical implementations, we have to use discrete solution Fd to 
approximate F. A complete justification of the numerical scheme requires us to 
prove \\(f) o F ~ (f) o FiiW —* as we refine the mesh. In [3], Allaire and Brizzi 
proved the convergence with respect to the discrete map F^ in the periodic case 
using asymptotic expansion, as well as some regularity assumptions requiring the 
mappings F and Fa smooth enough. However, in the general case, neither the tool 
of asymptotic expansion nor smoothness assumption is available, which makes the 
complete justification very difficult. Some discussions and further suggestions for 
similar problems in the context of variational mesh generation can be found in [17] . 
Another problem is, although F is guaranteed to be invertible, Fh is not. Fortu- 
nately this can be relieved if F is solved by piecewise linear finite element and the 
mesh only has non-obtuse-angled triangles |16j. In view of the above discussion, 
we need to compute Fa at a fine mesh such that the error ||F — F^^H is very small 
to fully resolve the small scale structure of F. 

We use the following notation 



a[v,w]:= / ^\/v{x,t)a{x)'Vw{x,t) dx. (2-11) 

Jn 

For V G Hq{Q) write TZh,oV the Ritz-Galerkin projection off on V^ with respect 
to the bilinear operator a[-,-], i.e., the unique element of V^ such that for all 

a[w,v-nhflv] = 0. (2.12) 

Define Y^ the subspace of L'^{0,T; H^{n)) as 

Y^ := {v G L^{0,T;H^{n)) : v{x,t) G V'^^t G [0,T]}. (2.13) 

Write Uh the solution in 1^ of the following system of ordinary differential 
equations: 

{K~^ip{x), dfuh)L'^n) + a['ip{x),Uh] = ii'ix),g)L2(n) for all t G (0, T) and ip E V^' 
Uh{x,0) = nh^ou{x,0), 

dtUh{x, 0) = TZhfidtu{x, 0). 

(2.14) 
The following theorem shows the error estimate of the semidiscrete solution. 
We need more smoothness on the forcing term g and the initial data than As- 
sumption [211] to guarantee the 0{h) convergence of the scheme 02.141) . On the 
other hand, we can see that even if g and all the initial data are smooth, with gen- 
eral conductivity matrix a(a;), we can merely expect u G L°°(0,T, iJ^(r2)) instead 
of the improved regularity L'^{0,T,H'^{Q)) in the harmonic coordinates, and the 
convergence rates will deteriote for the conventional finite elements. 



Assumption 2.2. Assume that the forcing term g satisfies d^g G L'^{Qt), dtQ ^ 
L°°(0,r,L2 ((])), initial value m(x, 0) and dtu{x,Q) satisfy d'fu\x,Q) E H^ln) and 
Va{x)Vdtulx,0) e L^{Vl) or equivalently dlu{x,Q) G /.^(fi). 

From now on we will always suppose without explicitly mentioning that As- 
sumption 12.21 is satisfied in the discussion of numerical homogenization method. 

Theorem 2.2. Suppose that Cordes condition \2.1\ and Assumption \2.^A hold, we 
have 

\\dt{u - Uh){.,T)\\^^^^^ + ||(« - "'^)(-'^)Li(f^) - <^^(ll^*^IU-(o,r,L2{n)) + WdlgW^Ur) 

+||a,\(x,o)||Hi(Q) + ||5fM(x,o)|U2(f,)). 

(2.15) 

The constant C depends on Cx, n, Q, fi„, Kmm, K^g^^, amin; Omax; and 

||(Trace[a])-i||^^(^). 

2.3 Numerical Homogenization in Time and Space 

Let M E N, (tn = njj)o<n<M is a discretization of [0, T]. (ipi) is a C^ basis of X^. 
Write trial space Z!^ the subspace of Y^ such that 

Z!^ = {w eY,^ : w{x,t) =y^Ci{t)ipi(^F{x)), Ci{t) are linear on (t„,t„+i] and 

i 

continuous on [0,T]} 

(2.16) 

Let test space U^ be the subspace of Y^ such that 

Ut = {iP EY^ : tpix.t) = ^diipi{F{x)), di are constant on [0,T].}. (2.17) 

i 

Write Vh the solution in Z^ of the following system of implicit weak formulation: 
forn e {0, . . . , M - 1} and V' e t/^^, 



ftn+l 

{K-'^lj,dtVh){tn+l) - {K~'^,dtVh){tn) = I {K-'dt^,dtVh)dt 



n + 1 



(2.18) 



a[ilj,Vh]dt+ / {'ip,g)dt. 



In equation (I2.18p . dtVhit) stands for lime|o(w/i(^) — Vh{t — e))/e. Once we 
know the values of Vh and dtVh at t„,, (12.181) is a linear system for the the unknown 



coefficients of dtVh(tn+i) in V'^. By continuity of f^ in time, we can obtain Vh{tn+i) 
by 

Vh(tn+l) = dtVh{tn+l){tn+l " ^n) +Vh{tn). (2.19) 

The following Theorem 12.31 shows the stability of the implicit scheme (12.181) : 

Theorem 2.3. Suppose that Cordes condition \ 2. 1\ and Assumption \2.2\ hold, we 
have 

\\dtVh{;T)\\L2^^) + \\vh{.,T)\\Hi^^) < C{\\g\\L2(^aT) + \\dMx,0)\\L2(^n) 

+ \\uix,0)\\Hi{n)) + Ch{\\dtg\\L'^(o,T,L2(n)) + Wd^gWL^nT) 
+ \\dM^,0)\\mm + ||9^(x,0)|U2(f,)). 

(2.20) 

The constant C depends on a^^^, a^aax, -K^min; -^max; o^nd T. 

The following Theorem 12.41 gives us the error estimate for the scheme (I2.18p . 

Theorem 2.4. Suppose that Cordes condition \2.1\ and Assumption \2.2\ hold, we 
have 

\\{dtUh-dtVh){.,T)\\L2^^n) + \\{uh-Vh){.,T)\\^,^^^<C^t{l + h~^) 

(II%I|l-(o,t,l2(c)) + WdhWh^^T) + ll'9i^w(a^>0)||^i(n) + ||(9fM(x,0)||L2(f^)). 

The constant C depends on Cx, T, amm; Omax? -^min; -^max; fJ'cr, o,nd 

||(Trace[a])-i||^^(^). 

3 Proofs 

The proofs are organized into three subsections corresponding to the three subsec- 
tions of section [21 

3.1 Compensation Phenomena: Proof of Theorem 12JJ 
Lemma 3.1. We have 

WdMhmiT) + a[dM{T) <CiT,^,K^^^)(a[d,um 

Amir, V /g j^N 



'-mm 



+ \\dM^.mhin) + \\dt9\\l2^nr 



Proof. In case a is smooth, differentiating (11. ip with respect to t, we have 

R-^dfu - div aVdtU = dtg. (3.2) 

multiplying by d^u, and integrating over Q, we obtain that 

2^ll^~^^t^«lli2(c) + 2^"['^*"] = (^*^' ^t^")i'(f2)- (3-3) 

Integrating the latter equation with respect to t and using Cauchy-Schwartz 
inequality we obtain that 



\K-^dMh^n)iT) + a[dM{T) <\\K-WMlHn)iO) + aM(0) 



(3.4) 
+ \\dt9\\LHnT)\\9tu\' 



Consider the following differential inequality, suppose that A is constant, 
B{t) > and nondecrease, X{t) > and X{t) is continuous with respect to 
t, 

X{t)<A + B{t){ X{s)dsY. (3.5) 

Jo 

Write Y{t) = sup^g^ ^j X{s), one has 

X(t) <A + B{t)t^ {Y{t)) ^ < A + ^iB{t)f + Y{t) ^ ^g g^ 

Take the supremum of both sides over t G [0, T], we have 

Y{T)<2A + T[B{T)Y. (3.7) 

It follows that 

K ■ 

+ ll5Mli^(n)(o) + ll%lli.(n, 



|5Mli^(n)m+«Mm<Cm§^,^max)(a[9,«](0) 



V) 



(3.8) 



In the case where a is nonsmooth we use Galerkin approximations of u in (II. ip 
and then pass to limit. This technique is standard and we refer to [IHl Section 
7.3. 2. cl for a reminder. D 



Lemma 3.2. 

l|5t«lli^(n)(2^)+«Mm<C^m^,^max)(aM(0) 



(3.9) 



Proof. Multiplying (11. ip by dtu, and integrating over Q, we obtain that 

2^ll^"^^t^lli2(f^) + 2 dt"[^] " ^^' '^*^)^'(f^)- (3.10) 

The remaining part of the proof is similar to the proof of Lemma 13.11 D 

We now need a variation of Campanato's result [H] on non-divergence form 
elliptic operators. For a symmetric matrix M, let us write 

Trace (M) 
"^^ - Trace(*MM)- ^^'^^^ 

Consider the following Dirichlet problem: 

Lmv = f (3.12) 

with Lm := J2ij=i Mij{x)didj. The following Theorem 13.11 is an adaptation of 
Theorem 1.2.1 of pi]. They are proved in [21] under the assumption that M is 
bounded and elliptic. It can be proved that the conditions fiM < oo and um < oo 
are sufficient for the validity of the theorem, we refer to [27] and [2H] for that proof. 

Theorem 3.1. Assume that ^m < oo, um G L°°{Q) andQ is convex. If f G L^(fi) 
the Dirichlet problem 1^3.12) has a unique solution satisfying 

11^1^2.2(0 ^ CfiMlWMfhHQ)- (3.13) 

Remark 3.1. The theorem can be extended to dimension n > 2 under the general 
Cordes condition [2] 



Let us now prove the compensation result in Theorem 12. 1[ Choose 

M:= ^oF-\ (3.14) 

|det(VF)|^ 

Recall that a := *VFaVF. (I3.14p is well defined since /xm = A^o- and 

\Wm\\1^(q) < T T-r II (Trace [a]) -^11^^ (3.15) 



Fix tG [0,T]. Choose 



|det(VF)|i 



(3.16) 



By the change of variable y = F{x), one obtains that 

|L2(n) < 2K^JdMLHu) + 2\\g\\L2^ny (3.17) 



10 



Using the notation k{y) := K{F ^{y)), g{y,t) := g{F ^iy),t), and u{y,t) := 
u{F~^{y),t), it follows from Theorem 13.11 that there exists a unique v G W'^''^{Q) 
such that 

J2{aiF-\y)))^fid,viy,t)=K-\y)dMy,t)-~9iy,t), (g^g) 

and 

||'y||w'2,2(f^) < CfiM\WM\\L°°{n){K^iJ\d^u\\L2(^Q) + ||5'||l2(c))- (3.19) 

By change of variable y = F{x) and the identity div aVF = we deduce that 
fl3.18p can be written as 

div {aV{v o F)) = K'^d'^u - g. (3.20) 

If d^u G L'^{Q) and g{.,t) G L'^{Q) we can use the uniqueness property for 
the solution of the following divergence form elliptic equation (with homogeneous 
Dirichlet boundary condition) 

div (aVw) = K-^d'fu - g. (3.21) 

to obtain that v o F = u. Thus we have proven Theorem 12.11 

3.2 Numerical Homogenization in Space: Proof of Theo- 
rem Eia 

In the following sections we will prove the convergence of semidiscrete and fully 
discrete numerical homogenization formulation 02.141) and 02.18p . 

We have the following lemmas which are the discrete analogs of Lemma 13.11 
and Lemma [3.41 

Lemma 3.3. We have 



II^Xlli2(n)(^)+«[5tW/.](r) <C(2^,T^,^max) a[9iM,](0) 



Lemma 3.4. 

\\dtu4l,^^){T) + a[uu]{T) <C(T,^,i^^ax)(aK](0) 

11 



(3.22) 



(3.23) 



Write TZh the projection operator mapping L^(0, T; HliVt)) onto Y^, such that 
for all V eY^: 

AT[v,u-nhu]=0 (3.24) 

let p := u — TZhU and 6 := T^/jit — m/^, where u/j is the solution of fl2.14p . 
For fixed t G [0, T] and v e ifQ(r2), we write 7lh,tv{-,t) the solution of: 

V^a(a;) V(t; - 7^?,,i^;(x, t)) da; = for all ij e V" (3.25) 

It is obvious that 7lhu{.,t) = 7lh,tu{.,t). For example, we can choose a series of 
test functions in (I3.24p which is separable in space and time, v{x,t) = T{t)X{x), 
T{t) is smooth in t and has S{t) function as its weak limit. 

We need the following lemma: 

Lemma 3.5. For v G Hq{Q) we have 

1 1 j_ 

{a[v -TZh,tv])^ ^ Cha^i,^n^\\v\\w2,2 (3.26) 

Proof. Using the change of coordinates y = F{x) we obtain that (write v := 
voF-^) 

a[v] = Q[d] (3.27) 

with 

Q[w] := / 'Vw{y)Q{y)Vw{y) dy (3.28) 

Jn 

and 

Using the definition of TZhV we derive that 

Q\v-nhVoF'H= inf Qlv-'A. (3.30) 

By interpolation property 02.71) it follows, 

Q[v - UhV o F-i] < A^ax(Q)Ci/i2||i)f^2,2(^). (3.31) 

where Amax(Q) is the supremum of eigenvalues of Q over Vt. 
It is easy to obtain that 

1 

Amax(vj — ^(^maxpcr [O.OZ) 

which finishes the proof. D 

12 



We will use the Lemmas 13.51 13.61 13.71 13.81 and 13.91 to obtain the approximation 
property of the projection operator TZh- 

With the improved Assumption 12.21 differentiate (11.11) with respect to t, and 
follow the proof of Theorem 12.11 we have 



Lemma 3.6. dt{u o p-^) G L°^{0,T,H^{n)) and 

\\dt{u o F'^)\\Lo.^o,T,H^n)) <C(||%||L-(o,r,L2(n)) + \\dt9\\mnT) 

+ ||9fu(x,0)|U2(n) + \\dMx,0)\\HHn)). 

The constant C is the one given in Theorem \2.1[ 

Apply Lemma [3.51 to dfU, we have 

Lemma 3.7. 

{AT[dtp])^ <Ch{\\dtg\\L--{o,T,L2{n)) + WdtgHmnr) 
+ \\dMx,0)\\LHn) + \\dMx,0)\\min))- 



(3.33) 



(3.34) 



The constant C depends on Cx, n, Q, fi^, ajam, ci^a.x, -^minj-^max; o,nd 

||(Trace[a])-i||^^(^) 

We have the following estimate for ||i9fp|| using the so-called Aubin-Nitsche 
trick [5]. . 

Lemma 3.8. 

Wdtph^inr) <Ch^{\\dtg\\L^io,T,L\n)) + Wd^gh^nr) ,^ ^^. 

+ ||9^(x,0)|U2(n) + ||9^(x,0)||Hvn)). 



The constant C in Lemma depends on Cx, n, Q, fi„, a^i^, flmax; -^min; -^max? O'l^'d 

lL°°(n) 



(Trace [o"]) "^1 



Proof. We choose v G L'^{0,T,HQ{fl)) to be the solution of the following linear 
problem: for all w G L^{0,T, H^{n)) 

At[w,v] = {w,dtp)L2(nj,). (3.36) 

Choosing w = dtp in equation (I3.36P we deduce that 

\\dtp\\l.^^^^)) = AT[dtp,v -Tluv]. (3.37) 

Using Cauchy Schwartz inequality we deduce that 

WdtPWh^nr) < {AT[dtp\Y{AT[v-nuv]Y. (3.38) 

13 



Since dtp{-,t) G L'^{ft), applying Theorem 13.11 for t G [0,T] then integrate over t, 
we obtain that 

\\v\\L2(o,T,w2^^n)) < C\\dp\\L2(nT)- (3-39) 

Using Lemma 13.51 we obtain that 

{At[v - 7^;,^;]) ^ < Ch\\dtp\\mnr)- (3-40) 

It follows that 

\\dtp\\LHnr)<Ch{AT[dtp])~'. (3.41) 

We deduce the lemma by applying Lemma [H751 to bound Axldtp]. D 

We have the following estimates for initial data, 
Lemma 3.9. 

\\nh,odtu{x,0) -dtu{x,0)\\L^n) < Ch'^{\\dtg{x,0)\\L2(^n) + \\d'fu{x,0)\\Hi{n) + \\dtu{x,0)\\L^n)) 

||7^ft,,ou(a;,0) - u{x,0)\\H^(^n) < Ch{\\dtg{x,0)\\L^n) + \\dfu{x,0)\\m{n) + \\dtu{x,0)\\L^n)) 

(3.42) 

Proof. We can estimate ||i9ip||L2(f^) using the duality argument similar to Lemma 
13.81 and derive the second inequality by Lemma [3.51 D 

Lemma 3.10. we have 

\\dt{u - M/i)||i2(f^)(T) + a[u - Uh]{T) < C{K^in,K^s,^,T)(\\dt{u - m/,)||^2(s^)(0) 

+a[u - Uh]{0) + \\dtp\\L^nT)\\dt{u - Uh)\\mnT) + ^r[<9tp] 

(3.43) 
Proof. For ^ G L'^{0,T,H^{n)), we have 

(i^-V, d^{u -Uh)) + a[ij, u-Uh] = 0. (3.44) 

Let ip = dfO = dtiu — Uh) — dtp, it follows 

2 -iTll^~^'9t(M - Uh) 11^2(0) + 2^^[^ ~ ^h] = (K'^dtp, df{u - Uh)) + a[dtp, u-Uh]- 

(3.45) 
Integrate with respect to t, using Cauchy-Schwartz inequality, we have 

^\\K--2dt{u - ix,)||i2(n)(T) - l\\K-'Wt{u - Uh) \\h^n)iO) + ^^u - w,](T) 



1 /-r _ 

- -a[M-M/,](0) < / Kj„||9tp||L2(n)||9j2(M-Mft)||i2(f^)rft+ [AT\dtp\AT\u - Ui^^ . 

(3.46) 
The remaining part of the proof is similar to the proof of Lemma 13.11 D 

Theorem 12.21 is a straightforward combination of Lemma 13. H Lemma 13. 3[ 
Lemma [3. 7t Lemma [3. 8^ Lemma [3.91 and Lemma [3.101 
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3.3 Numerical Homogenization in Space and Time: Proof 
of Theorems [23] and [274] 

Stability Choose ip E U^ in equation (12.181) such that ip{x,t) = dtVh{x,t) for 
t G {tn,tn+i]- We obtain that 

^ rtn+i 

\\K~^9tVh\\l2(^^){tn+i) - {K~^dtVh{tn+i) , dtVh{tn))L^n) = " / a[dtVh,Vh]dt 

rtn+i 

+ / {dtVh,g)mn)dt. 

(3.47) 
Observing that 

/ a[dtVh,Vh]dt = -a[vh\{tn+i) - 7^«K](^n)- (3.48) 

using Caucliy-Schwartz inequality it follows, 

ptn + l 

+2 / {dtVh,g)L^(n){t)dt. 
Summing over n from to M — 1, we have, 



(3.49) 



\\K-'^dtVh\WT) + a[vh]{T) < \\K-'^dtVh\\''{Q) + a[vhm + 2 f {dtVh,g)mn)dt. 

Jo 

(3.50) 

We conclude the proof of Theorem 12.31 using the inequality (13. 7p in the proof 

of Lemma 13. 1[ 

H^ Error Estimate We derive from equations (12.181) and (12.141) that 

{K'^ip, dtUh - dtVh){tn+i) - (K'^ip, dtUh - dtVh){tn) 
ptn+i ptn+i (3,51) 

- / [K dtip, dtUh - dtVhjdt + / a[ijj, Uh - Vh]dt = 0. 

Let ip = dfith — dfVh where Uh is the linear interpolation of Uh over Z!^. Write 
Vh = Uh — Vh and Wh = Uh — Uh, it follows that 

{K-'^dtyh,dtyh){tn+i) + {K-^dtWh,dtyh){tn+i) - {K-^dtyh,dtyh){tn) 

rtn + l /-tn + l (3,52) 

-{K dtWh,dtyh){tn) + a[dtyh,yh]dt+ a[dtWh,yh]dt = 0. 



in " i; 
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Observing J^ """"^ dtWh{x, t) dt = we need the following lemma, which is a slight 
variation of the Hilbert- Bramble lemma, |10j . 

Lemma 3.11. If J^^^' u{s) ds = 0, then 

u^ < jAt I ^ u'isfds. (3.53) 

Since dfwh{x,t) = —d^Uh{x,t) in (t„,t„+i], by Lemma 13.111 we have 

/ \dtWh{x,t)\^dxdt<-At / / \dtUh{x,t)\^dxdt, (3.54) 

and 

/ \dtWh{x,t)\^dxdt<-At^ \d'fuh{x,t)\'^dxdt. (3.55) 

Jt„ Jn 4 Jf^ J^ 

Using the inverse Sobolev inequality (12. 9p we obtain from equation (I3.55P that 

/ / a[dtWh]dxdt < C-— / / \dfuh{x,t)\'^ dxdt. (3.56) 

Jt„ Jn "■ Jt„ Jn 

Summing (I3.52p over n, notice yh{0) = 0, dtyuiS^) = we obtain that 

1 fT 

{K ^dtyh,dtyh)L2(n){T) + -a[yh{.,T)] = - a[dtWh,yh]dt - {K \dtWh,dtyh)LHn){T). 

(3.57) 

Theorem 12.41 is a straightforward consequence of (13.570 , the estimates (I3.54p , 
(I3.56p . Lemma EH] and Lemma EiS 

4 Numerical Experiments 

In this section, we will present the numerical algorithm and examples. 

We use web extended B-spline based finite element [IH] to span the space X'^ 
introduced in subsection 12. 2[ For all the numerical examples, we compute the 
solutions up to time T = 1. The initial condition is u{x,0) = and Ut{x,0) = 0. 
The boundary condition is u{x,t) = 0, x G dQ. For simplicity, the computational 
domain is the square [—1, 1] x [—1, 1] in dimension two. 

We have a fine mesh and a coarse mesh characterized by different degrees of 
freedom (dof). In general, the fine mesh is generated by hierarchical refinement of 
the coarse mesh: for each triangle of the coarse mesh, choose middle points of its 
3 edges as new vertices, and divide the triangle into 4 new triangles, a is defined 
as a piecewise constant function over each fine mesh triangle, and is evaluated at 
the center of mass of the triangle. 
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Algorithm 4.1 (Algorithm for Numerical Homogenization) . 

1. Compute F on fine mesh, the fine mesh solver for F is Matlah routine assempde. 

2. Construct multi-scale finite element basis ip = if o F, compute stiffness matrix K and 
mass matrix M. 

3. March (12.181) and (12.191) in time with respect to the coarse dof. 

4. Repeat 3 if we have multiple right hand sides. 

In the implementation, F is approximated by a piecewise linear finite element 
solution. We mesh the square such that no triangle has an obtuse angle, therefore 
F is an invertible piecewise linear mapping [16]. When we construct tp, we simply 
take its piecewise linear interpolation on the fine mesh. 

All the computations were done at a single Opteron Dual-Core 2600 cpu of a 
Sun Fire X4600 server, and programmed in Matlab 7.3. 

Example 4.1. Multiscale trigonometric coefficients 

The following example is extracted from [2S] as a problem without scale sepa- 
ration: 



1 /l.l + sin(27rx/ei) 1.1 + sin(27ry/e2) 1.1 + cos(27rx/e3) 



6 V 1.1 + sm{2ny/ei) 1.1 + cos(27ra;/e2) 1.1 + sin(27r?//e3) 

1.1 + sin(27ry/e4) 1.1 + cos(27ra;/e5) . .. 2 2n , . 

H ; ; — r H ; ; — r + sm 4x v ) + I 

1.1 + cos(27ra;/e4) 1.1 + sm{2ny/e5) ^ ^ ' 



(4.1) 



where ei = |, €2 = ^, €3 = j^, £4 = ■^, £5 = ^. The conductivity a is smooth, 
therefore it satisfies Cordes condition 12.11 

First, we want to compare the performance of different numerical homogeniza- 
tion methods, 

• LFEM: A multi-scale finite element where F is computed locally (instead 
of globally) on each triangle K of the coarse mesh as the solution of a cell 
problem with boundary condition F{x) = a; on OK. This method has been 
implemented in order to understand the effect of the removal of global infor- 
mation in the structure of the metric induced by F. 

• Y'EM.jipiin. The Galerkin scheme using the finite elements ipi = cpioF, where 
ifi are the piecewise linear nodal basis elements. 

• FEMjipsp'- The Galerkin scheme using the finite element ipi = ipiO F, where 
ipi are weighted cubic B-spline elements. 

Suppose Uf is the finite element solution of (ll.lj) computed on the fine mesh 
at time T = 1, the fine mesh solver is Matlab routine hyperbolic, which uses linear 
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Table 1: Example 14.11 numerical errors of different methods, coarse dof 49, fine 
dof 261121, g = 



Method 



LFEM 

FEM_lljiin 

FEM_ijsp 



L' 



0.0440 
0.0315 
0.0021 



0.0982 
0.0518 
0.0035 



L' 



0.0534 
0.0362 
0.0022 



H' 



0.2054 
0.1601 
0.0189 



Table 2: Example 14.11 numerical errors of FEM^ipsp, with g = 1, doff is fine mesh 
dof, dofc is coarse mesh dof 



doff 


dofc 


Li 


^oo 


L' 


H' 


65025 


9 

49 
225 


0.0075 
0.0023 
0.0009 


0.0118 
0.0037 
0.0023 


0.0074 
0.0023 
0.0010 


0.0394 
0.0194 
0.0117 


261121 


9 

49 
225 


0.0070 
0.0021 
0.0009 


0.0106 
0.0035 
0.0025 


0.0069 
0.0022 
0.0010 


0.0373 
0.0188 
0.0117 



finite element basis in space and adaptive ODE integrator in time, v^ is the 
solution of (12.181) . Numerical errors in the norm || • || are computed by 



error 



\Vh 



Uf 



(4.2) 



Numerical errors in L^, L^, L°° and H^ norm are computed. 

In Table [H performances of different methods with coarse mesh dof 49 are 
compared. We observe that the methods using global F have better performance, 
and FEM_ipsp is much better than other methods. Note that the improvement of 
FEM jipiin over LFEM is not as significant as the elliptic case [28]. 

From now on, all the results are computed by the method FEMjip^p. 

Next, the impact of right hand side on accuracy will be investigated. We solve 
equation (11.11) with a time independent source term g = 1, a. slowly varying term 
g = sin(2.4x — 1.8y + 27rt), and a Gaussian source term given by 



9{x,y) 



1 



V2 



TTcr^ 



a;2 + (^_ 0.15)2 



(4.3) 



with a = 0.05. Notice that as cr ^ 0, the source function will become singular in 
space. 
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Table 3: Example 14. 11 numerical errors oiFEM^ipsp, with g = sin(2.4a; — 1.8y+27rt) 



doff 


dofc 


L' 


L°° 


L' 


H' 


65025 


9 

49 
225 


0.0400 
0.0107 
0.0035 


0.0390 
0.0105 
0.0047 


0.0360 
0.0096 
0.0033 


0.0869 
0.0393 
0.0233 


261121 


9 

49 
225 


0.0399 
0.0104 
0.0034 


0.0373 
0.0109 
0.0047 


0.0359 
0.0095 
0.0033 


0.0866 
0.0391 
0.0231 



Table 4: Example 14.11 numerical errors of FEMjipsp, with the Gaussian source g 



m 



doff 


dofc 


Li 


L°° 


L' 


H' 


65025 


9 
49 

225 


0.0581 
0.0272 
0.0096 


0.2270 
0.1023 
0.0179 


0.0704 
0.0333 
0.0095 


0.3484 
0.2305 
0.0957 


261121 


9 

49 
225 


0.0574 
0.0274 
0.0097 


0.2199 
0.976 
0.0212 


0.0688 
0.0332 
0.0101 


0.3436 

0.2254 
0.1005 



Table [2] presents the errors for time independent term g = 1. Table [3] presents 
the errors for slowly varying source term g = sin(2.4x — 1.8y + 27rt). Table H] 
presents the errors for relatively singular Gaussian forcing term. It is clear that 
if the source term is time independent and smooth in space, the method is more 
accurate, which corresponds to the smoothness requirement of g in Theorem 12.41 
In all the examples, we have tried two fine meshes which have dof 65025 and 
261121 respectively, roughly 250 x 250 and 500 x 500. It can be seen that for fixed 
coarse dof, the errors with respect to different fine dof are pretty close, which 
means fine dof 65025 is enough for this problem. 

The Figure [1] shows the L^ error evolution with respect to time, which is typ- 
ical for other norms. The overshoot at the beginning corresponds to the time 
discretization step. After several steps, the errors tend to be stable. 
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Figure 1: Example l4.lt L^ error with respect to time, for g = 1, coarse dof 9, fine 
dof 65025 



Example 4.2. Time independent high conductivity channel 

High conductivity channel is an interesting test problem in many petroleum ap- 
plications because of its strong non-local effects. In this example a is characterized 
by a narrow and long ranged high conductivity channel. We choose a{x) = A^ 1, 
if X is in the channel, and a (x) = 1, if a; is not in the channel. The media is illus- 
trated in Figure [21 However, in this case, whether or not Cordes condition 12.11 is 
not clear. We will go ahead testing the numerical performance of our method. 




Figure 2: Example 14.21 high conductivity channel medium 



Table [3] shows numerical errors for g = 1 with fixed coarse dof 49 and A = 
10^, 10^, 10'^, 10^ respectively. From the table we can see that the errors grow with 
the aspect ratio increasing, but the growth is quite moderate and the numerical 
behavior of the method is stable. The errors for time dependent right hand side 
g = sin(2.4x — 1.8y + 27rt) with A = 10^ are also given in Table El 
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Table 5: Example 14.21 numerical errors with respect to different aspect ratios, 



dof 261 : 


L21. 








A 


L' 


L°° 


L^ 


H' 


10 


0.0021 


0.0056 


0.0025 


0.0240 


100 


0.0118 


0.0497 


0.0180 


0.0964 


1000 


0.0181 


0.0931 


0.0316 


0.1308 


10000 


0.0243 


0.1174 


0.0419 


0.1550 



Table 6: Example 



4.21 numerical errors for g = sin(2.4a; — 1.8y + 27rt) 



doff 



65025 



261121 



dofc 



9 

49 
225 



9 

49 
225 



0.0750 
0.0301 
0.0096 



0.0752 
0.0302 
0.0094 



0.0777 
0.0351 
0.0118 



0.0779 
0.0345 
0.0116 



0.0729 
0.00298 
0.0092 



0.0731 
0.0299 
0.0091 



H' 



0.1528 
0.0779 
0.0324 



0.1533 
0.0782 
0.0321 



Example 4.3. Time independent site percolation 

In this example we consider the site percolating medium associated to Figure [3l 
In this case, we subdivide the square into a 64 x 64 checkerboard, the conductivity 
of each site is equal to 7 or I/7 with probability 1/2. We have chosen 7 = 10 in 
this example. In fact, this medium may not satisfy the Cordes condition 12.11 (also 
refer to Remark 12. ip . However, we will show that the method still works fine for 
this example. 

Figure m shows u computed with 261121 dof and Vh computed with 9 dof in 
the case g = 1 at time 1 using method FEM jipgp. They are visually almost the 
same even in terms of small scale features. Table [7| gives the numerical errors for 
g = I with respect to different coarse and fine dof 

Finally, we consider the site percolating medium, with Neumann boundary 
condition and a more realistic forcing term. The source term is given by g{x, t) = 
T{t)X{x,y), X{x,y) is the Gaussian source function described by 



X{x,y) = 
with a = 0.05, T{t) = Ti(t)T2(t) 



V2 



TTO"^ 



:exp ( 



x'^ + y'^ 
2^2 



)> 



10 



Ti(t)=^2 



1- 



kn 



sin{2kTTt), 



(4.4) 



(4.5) 
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Figure 3: Example 14. 3[ site percolation medium 



Table 7: Example 14. 3^ numerical errors for g = 1 



doff 


dofc 


L' 


L~ 


L^ 


H' 


65025 


9 

49 
225 


0.0750 
0.0301 
0.0135 


0.0777 
0.0351 
0.0147 


0.0729 
0.00298 
0.0133 


0.1528 
0.0779 
0.0333 


261121 


9 

49 
225 


0.0752 
0.0302 
0.0131 


0.0779 
0.0345 
0.0145 


0.0731 
0.0299 
0.0130 


0.1533 
0.0782 
0.0329 



and T2(t) = erfc(8(t — 0.5)), erfc is the complementary error function. We use 
this source term to emulate a source acting around the origin before t = 0.5, then 
suddenly decays. See Figure Ofor T(t) in (0, 1). 

In fact, our future goal is to simulate the response of an explosion, usually this 
is done with a so-called Ricker function [6] , i.e., g{x, y, t) = 6o{x — Xs,y — ys)R(t) 
with 

Rit) = (1 - 2n\fot - 1)2) exp(-7r2(/ot - 1)^] (4.6) 

where ^o is the Dirac function and /o is called the central frequency of the source 
wavelet. It is clear that Ricker function does not belong to L'^{Qt), our analysis 
does not apply and numerical experiment failed in this case. Therefore we would 
like to test the above modified source term first. 

Numerical errors for this modified source are given in Table [HI The errors 
are acceptable but not so good as g = 1. Possible solutions include the adaptive 
integration in time and adaptation in space around the source. 
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(a) 



(b) Vh. 



Figure 4: u computed with dof 261121 and Vh computed with dof 9 at time 1, 
they are interpolated on a coarser mesh in order to have a clear picture. 




Figure 5: t — > 5'(0,t) 



Table 8: Example 14.31 numerical errors for modified source 



doff 


dofc 


L^ 


^oo 


U' 


H^ 


65025 


9 

49 
225 


0.0429 
0.0261 
0.0185 


0.1580 
0.0965 
0.0879 


0.0505 
0.0324 
0.0226 


0.4065 
0.3030 
0.2718 


261121 


9 

49 
225 


0.0491 
0.0263 
0.0178 


0.1189 
0.0969 
0.1139 


0.0578 
0.0324 
0.0221 


0.4259 
0.2780 
0.2474 
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5 Conclusion and Further Remarks 

From above analysis and numerical examples, we observe that good numerical 
approximations can be obtained with much fewer degrees of freedom for acoustic 
wave equation with heterogeneous coefficients, even for the cases which do not 
satisfy the Cordes condition. Compared with the multiscale finite element method 
which compute the basis locally, our method has much better accuracy, especially 
for problems with strong non-local effects. 

As it has been done in [2H], once one understand that the key idea for the 
homogenization of f 1 1.11) lies in its higher regularity properties with respect to har- 
monic coordinates one can homogenize (11. ip through a different numerical method 
(such as a finite volume method). 

Moreover, it could be observed that one could use any set of n linearly inde- 
pendent solutions of (11.11) instead of the harmonic coordinates. The key property 
allowing the homogenization of (11.11) lies in the fact that if the data (right hand 
side and initial values) has enough integrability then the space of solutions is at 
small scales close in H^ norm to a space of dimension n. Thus once one has ob- 
served at least n linearly independent solutions of (II. ip , one has seen all of them 
at small scales. 

Write L := -VaV. L'^ maps H-^{Q) into H^{Q), it also maps L'^{Q) into V 
a sub- vector space of Hq{Q). The elements of V is close in H^ norm to a space of 
dimension n (the dimension of the physical space Q) in the following sense. 

Let 7/i be a triangulation of ^2 C MJ^ of resolution h (where < h < diam(i7)). 
Let A set of mappings from T^ into the unit sphere of M"'^^ (if A G A then Ais 
constant on each triangle K eT^ and ||A(7^)|| = 1), then 

II v-^"-+l \ II 

sup ml --j < Ch 5.1 

vi,V2,...,v^+i(iV ^(^^ l^i^l \\ya\/Vi\\L2^p_) 

Equation (15. ip is saying that any n+1 elements of V are (at an h approximation 
in H^ norm) linearly dependent. Recall that n + 1 vectors are linearly dependent 
in a linear combination (with non null coefficients) of these vectors in the null 
vector. In (15.11) the linear combination of the n + 1 vectors is at relative distance 
of order h (resolution of the triangulation) from 0. 

We notice that some recent results using global information [U [22l [21] are 
formulated in a partition of unity framework [^. In this case, {1, Fi, ■ ■ ■ , -F„} can 
be used to construct the local approximation space. 
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